function dy = scaffold_dynamics(t,y,P)

P = num2cell(P);

[k1f,k1r,k2,k3f,k3r,k4,k5f,k5r,k6,k7f,k7r,k8f,k8r,k9f,k9r,k10f,k10r,k11] = deal(P{:});

y = num2cell(y);

[S,A0,A,B0,B,X,Y,SA0,XA,YB,AC,A0C,BC,B0C,ABC,A0BC,AB0C,A0B0C] = deal(y{:});

% Phosphorylation of A:
r1 = k1f*S*A0 - k1r*SA0;
r2 = k2*SA0;
% Dephosphorylation of A:
r3 = k3f*X*A - k3r*XA;
r4 = k4*XA;
% Dephosphorylation of B:
r5 = k5f*Y*B - k5r*YB;
r6 = k6*YB;
% Binding of C to A0 and A:
r7 = k7f*A*C - k7r*AC;
r8 = k10f*A0*C - k10r*A0C;
% Binding of C to B0 and B:
r9 = k8f*B0*C - k8r*B0C;
r10 = k9f*B*C - k9r*BC;
% Binding of B0C and BC to A0 and A:
r11 = k7f*A*BC - k7r*ABC;
r12 = k7f*A*B0C - k7r*AB0C;
r13 = k10f*A0*BC - k10r*A0BC;
r14 = k10f*A0*B0C - k10r*A0B0C;
% Binding of A0C and AC to B0 and B:
r15 = k8f*B0*AC - k8r*AB0C;
r16 = k8f*B0*A0C - k8r*A0B0C;
r17 = k9f*B*AC - k9r*ABC;
r18 = k9f*B*A0C - k9r*A0BC;
% Phosphorylation of B0:
r19 = k11*AB0C;

dS = -r1 + r2;
dA0 = -r1 + r4 - r8 -r13 - r14;
dA = r2 - r3 - r7 - r11 - r12;
dB0 = r6 - r9 - r15 - r16;
dB = -r5  - r10 - r17 - r18;
dX = -r3 + r4;
dY = -r5 + r6;
dSA0 = r1 - r2;
dXA = r3 - r4;
dYB = r5 - r6;
dAC = r7 - r15 - r17;
dA0C = r8 - r16 - r18;
dBC = r10 - r11 - r13;
dB0C = r9 - r12 - r14;
dABC = r11 + r17 + r19;
dA0BC = r10 + r18;
dAB0C = r12 + r15 - r19;
dA0B0C = r14 + r16;

dy = [dS;
    dA0;
    dA;
    dB0;
    dB;
    dX;
    dY;
    dSA0;
    dXA;
    dYB;
    dAC;
    dA0C;
    dBC;
    dB0C;
    dABC;
    dA0BC;
    dAB0C;
    dA0B0C];